

%% drop top 5% '09-'12 spending increases %%

dChar = readtable('../../../raw_data/district_characteristics.csv');
e09 = readtable('../../../raw_data/election_results_2009.csv');
e12 = readtable('../../../raw_data/election_results_2012.csv');
cs09 = readtable('../../../raw_data/campaign_spending_2009.csv');
cs12 = readtable('../../../raw_data/campaign_spending_2012.csv');

% order observations by candidate-coalition choice
d0 = find(e12.candCM==0); N0 = size(d0,1);
d1 = find(e12.candCM==1); N1 = size(d1,1);
d2 = find(e12.candCM==2); N2 = size(d2,1);
dChar = dChar([d0;d1;d2],:);
e09 = e09([d0;d1;d2],:); e12 = e12([d0;d1;d2],:);
cs09 = cs09([d0;d1;d2],:); cs12 = cs12([d0;d1;d2],:);
clear d0 d1 d2

% find largest spending discrepancies
sd = (cs12.PRI + cs12.PVEM + cs12.CM - ...
    (cs09.PRI + cs09.PVEM + cs09.PM)) ./ (cs09.PRI + cs09.PVEM + cs09.PM) <= ...
    prctile((cs12.PRI + cs12.PVEM + cs12.CM - ...
    (cs09.PRI + cs09.PVEM + cs09.PM)) ./ (cs09.PRI + cs09.PVEM + cs09.PM),95);
sd_long = [sd;sd;sd(1:N0);sd(N0+N1+(1:N2));sd(1:N0);sd(N0+(1:N1));sd];
cs12 = cs12(sd,:);
dChar = dChar(sd,:);
e09 = e09(sd,:);
e12 = e12(sd,:);
N0 = sum(e12.candCM==0);
N1 = sum(e12.candCM==1);
N2 = sum(e12.candCM==2); N = N0 + N1 + N2;

% district indicator/aggregation matrix
A = [1:N,1:N,1:N0,N0+N1+(1:N2),1:N0+N1,1:N]';
A = 1 * (A==A');

% lagged vote shares
P = [e09.MP, ...
    e09.NA, ...
    e09.PVEM + 0.5 * e09.PM, ...
    e09.PRI + 0.5 * e09.PM, ...
    e09.PAN] ./ e09.registered;

% 2012 vote shares
S = [e12.MP ./ e12.registered; ...
    e12.NA ./ e12.registered; ...
    e12.PVEM(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+N1+(1:N2),1) + e12.PVEM(N0+N1+(1:N2),1) + e12.CM(N0+N1+(1:N2),1)) ./ e12.registered(N0+N1+(1:N2),1); ...
    e12.PRI(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+(1:N1),1) + e12.PVEM(N0+(1:N1),1) + e12.CM(N0+(1:N1),1)) ./ e12.registered(N0+(1:N1),1); ...
    e12.PAN ./ e12.registered];


%% baseline specification

% demographics
D = [dChar.femaleHH,dChar.over60,dChar.rural];
dnames = {'female','over60','rural'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);
NN = size(X,1);

% add log-lagged vote shares
PP = zeros(NN,3);
PP([1:N0,N+(1:N0),2*N+(1:N0),2*N+N0+N2+(1:N0),3*N+N0+(1:N0)],1) = reshape(log(P(1:N0,:)),5*N0,1);
PP([N0+(1:N1),N+N0+(1:N1),2*N+N0+N2+N0+(1:N1),3*N+N0+N0+(1:N1)],2) = reshape(log([P(N0+(1:N1),1:2),sum(P(N0+(1:N1),3:4),2),P(N0+(1:N1),5)]),4*N1,1);
PP([N0+N1+(1:N2),N+N0+N1+(1:N2),2*N+N0+(1:N2),3*N+N0+N0+N1+(1:N2)],3) = reshape(log([P(N0+N1+(1:N2),1:2),sum(P(N0+N1+(1:N2),3:4),2),P(N0+N1+(1:N2),5)]),4*N2,1);
X = [X,sum(PP,2)];

% (endogenous) campaign spending
C = [cs12.MP;...
    cs12.NA;...
    cs12.PVEM(1:N0,1);cs12.CM(N0+N1+(1:N2),1);...
    cs12.PRI(1:N0,1);cs12.CM(N0+(1:N1),1);...
    cs12.PAN];
XC = [X,C,C.^2];
K = size(XC,2);


%% baseline specification: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S = [X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))]];

% OLS
Y2S = [log(e12.PVEM(N0+1:end)) - log(e12.CM(N0+1:end));...
    log(e12.PRI(N0+1:end)) - log(e12.CM(N0+1:end))];
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ ...
    (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table D2 (column (III)):')
disp(print_results(round(rB2S,3),dnames,[],0,[],1))


%% baseline specification: BLP

load('../../voting_stage.mat','xi0','rBblp','Z0')
xi0 = xi0(sd_long,:);
Z0 = Z0(sd_long,:);
KZ = size(Z0,2);

% SGI nodes and weights
[nu,wnu] = nwspgr('KPN',2,6);

% BLP
W = Z0' * diag(xi0.*xi0) * Z0; % GMM weighting matrix
Bblp = rBblp(:,1);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit5] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0;Z0'*xi0],[],[],...
    [sparse(KZ,length(Bblp)),Z0',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0' * Ghat;
VBblp = inv(Ghat'*((Z0'*diag(xi.*xi)*Z0)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp5 = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (III)):')
disp(print_results(round(rBblp5,3),dnames,[],0,1,0))


%% region FE

% demographics
D = [dChar.region<=2,dChar.region==3,dChar.femaleHH,dChar.over60,dChar.rural];
rfe = {'north','southeast'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);

% add log-lagged vote shares
X = [X,sum(PP,2)];

% (endogenous) campaign spending
XC = [X,C,C.^2];
K = size(XC,2);


%% region FE: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S = [X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))]];

% OLS
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S_RFE5 = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table D2 (column (VII)):')
disp(print_results(round(rB2S_RFE5,3),dnames,rfe,0,[],1))


%% region FE: BLP

load('../../voting_stage.mat','xi0rfe','rBblp_RFE','Z0rfe')
xi0rfe = xi0rfe(sd_long,:);
Z0rfe = Z0rfe(sd_long,:);
KZ = size(Z0rfe,2);

% BLP
W = Z0rfe' * diag(xi0rfe.*xi0rfe) * Z0rfe; % GMM weighting matrix
Bblp = rBblp_RFE(:,1);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit_RFE5] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0rfe;Z0rfe'*xi0rfe],[],[],...
    [sparse(KZ,length(Bblp)),Z0rfe',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0rfe' * Ghat;
VBblp = inv(Ghat'*((Z0rfe'*diag(xi.*xi)*Z0rfe)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_RFE5 = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (VIII)):')
disp(print_results(round(rBblp_RFE5,3),dnames,rfe,0,1,0))


%% drop top 10% '09-'12 spending increases %%

dChar = readtable('../../../raw_data/district_characteristics.csv');
e09 = readtable('../../../raw_data/election_results_2009.csv');
e12 = readtable('../../../raw_data/election_results_2012.csv');
cs09 = readtable('../../../raw_data/campaign_spending_2009.csv');
cs12 = readtable('../../../raw_data/campaign_spending_2012.csv');

% order observations by candidate-coalition choice
d0 = find(e12.candCM==0); N0 = size(d0,1);
d1 = find(e12.candCM==1); N1 = size(d1,1);
d2 = find(e12.candCM==2); N2 = size(d2,1);
dChar = dChar([d0;d1;d2],:);
e09 = e09([d0;d1;d2],:); e12 = e12([d0;d1;d2],:);
cs09 = cs09([d0;d1;d2],:); cs12 = cs12([d0;d1;d2],:);
clear d0 d1 d2

% find largest spending discrepancies
sd = (cs12.PRI + cs12.PVEM + cs12.CM - ...
    (cs09.PRI + cs09.PVEM + cs09.PM)) ./ (cs09.PRI + cs09.PVEM + cs09.PM) <= ...
    prctile((cs12.PRI + cs12.PVEM + cs12.CM - ...
    (cs09.PRI + cs09.PVEM + cs09.PM)) ./ (cs09.PRI + cs09.PVEM + cs09.PM),90);
sd_long = [sd;sd;sd(1:N0);sd(N0+N1+(1:N2));sd(1:N0);sd(N0+(1:N1));sd];
cs12 = cs12(sd,:);
dChar = dChar(sd,:);
e09 = e09(sd,:);
e12 = e12(sd,:);
N0 = sum(e12.candCM==0);
N1 = sum(e12.candCM==1);
N2 = sum(e12.candCM==2); N = N0 + N1 + N2;

% district indicator/aggregation matrix
A = [1:N,1:N,1:N0,N0+N1+(1:N2),1:N0+N1,1:N]';
A = 1 * (A==A');

% lagged vote shares
P = [e09.MP, ...
    e09.NA, ...
    e09.PVEM + 0.5 * e09.PM, ...
    e09.PRI + 0.5 * e09.PM, ...
    e09.PAN] ./ e09.registered;

% 2012 vote shares
S = [e12.MP ./ e12.registered; ...
    e12.NA ./ e12.registered; ...
    e12.PVEM(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+N1+(1:N2),1) + e12.PVEM(N0+N1+(1:N2),1) + e12.CM(N0+N1+(1:N2),1)) ./ e12.registered(N0+N1+(1:N2),1); ...
    e12.PRI(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+(1:N1),1) + e12.PVEM(N0+(1:N1),1) + e12.CM(N0+(1:N1),1)) ./ e12.registered(N0+(1:N1),1); ...
    e12.PAN ./ e12.registered];


%% baseline specification

% demographics
D = [dChar.femaleHH,dChar.over60,dChar.rural];
dnames = {'female','over60','rural'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);
NN = size(X,1);

% add log-lagged vote shares
PP = zeros(NN,3);
PP([1:N0,N+(1:N0),2*N+(1:N0),2*N+N0+N2+(1:N0),3*N+N0+(1:N0)],1) = reshape(log(P(1:N0,:)),5*N0,1);
PP([N0+(1:N1),N+N0+(1:N1),2*N+N0+N2+N0+(1:N1),3*N+N0+N0+(1:N1)],2) = reshape(log([P(N0+(1:N1),1:2),sum(P(N0+(1:N1),3:4),2),P(N0+(1:N1),5)]),4*N1,1);
PP([N0+N1+(1:N2),N+N0+N1+(1:N2),2*N+N0+(1:N2),3*N+N0+N0+N1+(1:N2)],3) = reshape(log([P(N0+N1+(1:N2),1:2),sum(P(N0+N1+(1:N2),3:4),2),P(N0+N1+(1:N2),5)]),4*N2,1);
X = [X,sum(PP,2)];

% (endogenous) campaign spending
C = [cs12.MP;...
    cs12.NA;...
    cs12.PVEM(1:N0,1);cs12.CM(N0+N1+(1:N2),1);...
    cs12.PRI(1:N0,1);cs12.CM(N0+(1:N1),1);...
    cs12.PAN];
XC = [X,C,C.^2];
K = size(XC,2);


%% baseline specification: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S = [X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))]];

% OLS
Y2S = [log(e12.PVEM(N0+1:end)) - log(e12.CM(N0+1:end));...
    log(e12.PRI(N0+1:end)) - log(e12.CM(N0+1:end))];
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ ...
    (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S10 = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table D2 (column (IV)):')
disp(print_results(round(rB2S10,3),dnames,[],0,[],1))


%% baseline specification: BLP

load('../../voting_stage.mat','xi0','rBblp','Z0')
xi0 = xi0(sd_long,:);
Z0 = Z0(sd_long,:);
KZ = size(Z0,2);

% SGI nodes and weights
[nu,wnu] = nwspgr('KPN',2,6);

% BLP
W = Z0' * diag(xi0.*xi0) * Z0; % GMM weighting matrix
Bblp = rBblp(:,1);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit10] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0;Z0'*xi0],[],[],...
    [sparse(KZ,length(Bblp)),Z0',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0' * Ghat;
VBblp = inv(Ghat'*((Z0'*diag(xi.*xi)*Z0)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp10 = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (IV)):')
disp(print_results(round(rBblp10,3),dnames,[],0,1,0))


%% region FE

% demographics
D = [dChar.region<=2,dChar.region==3,dChar.femaleHH,dChar.over60,dChar.rural];
rfe = {'north','southeast'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);

% add log-lagged vote shares
X = [X,sum(PP,2)];

% (endogenous) campaign spending
XC = [X,C,C.^2];
K = size(XC,2);


%% region FE: voting second stage

X2S = [blkdiag(ones(N1,1),ones(N2,1)),D(N0+1:end,:)];
X2S = blkdiag(X2S,X2S);
X2S = [X2S,[log(P(N0+1:end,3));log(P(N0+1:end,4))]];

% OLS
B2S = X2S \ Y2S;
xi2S = Y2S - X2S * B2S;
VB2S = (X2S' * X2S) \ (X2S' * diag(xi2S.*xi2S) * X2S) / (X2S' * X2S); % heteroskedasticity-robust variance
seB2S = sqrt(diag(VB2S)); % standard errors
rB2S_RFE10 = [B2S,seB2S,2*(1-normcdf(abs(B2S./seB2S)))];

disp(' ')
disp('for Table D2 (column (VIII)):')
disp(print_results(round(rB2S_RFE10,3),dnames,rfe,0,[],1))


%% region FE: BLP

load('../../voting_stage.mat','xi0rfe','rBblp_RFE','Z0rfe')
xi0rfe = xi0rfe(sd_long,:);
Z0rfe = Z0rfe(sd_long,:);
KZ = size(Z0rfe,2);

% BLP
W = Z0rfe' * diag(xi0rfe.*xi0rfe) * Z0rfe; % GMM weighting matrix
Bblp = rBblp_RFE(:,1);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit_RFE10] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0rfe;Z0rfe'*xi0rfe],[],[],...
    [sparse(KZ,length(Bblp)),Z0rfe',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0rfe' * Ghat;
VBblp = inv(Ghat'*((Z0rfe'*diag(xi.*xi)*Z0rfe)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_RFE10 = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (IX)):')
disp(print_results(round(rBblp_RFE10,3),dnames,rfe,0,1,0))


%%

save('voting_stage_dropOutliers.mat')



